LIBRARY  OF  THE 

UNIVERSITY  OF  ILLINOIS 

AT  URBANA-CHAMPAIGN 


510 /c4 
Cop.  2. 


Digitized  by  the  Internet  Archive 
in  2013 


http://archive.org/details/divergenceofston478diam 


it/* 

JL  Report  No.    478 


/ 


THE  DIVERGENCE   OF  STONE'S   FACTORIZATIONS 
WHEN  NO  PARAMETERS  ARE  USED 


by 


MARTIN  DIAMOND 


September,  1971 


1EUIBRARYOEIHI 


NOV     9  1972 


UNIVERSITY  OF  ILLINOIS 
AT  URBANA-CHAMPAIGN 


DEPARTMENT  OF  COMPUTER  SCIENCE 
UNIVERSITY  OF  ILLINOIS  AT  URBANA-CHAMPAIGN 


URBANA,   ILLINOIS 


Report  No.    1+78 


THE   DIVERGENCE   OF  STONE'S   FACTORIZATIONS 
WHEN  NO  PARAMETERS  ARE  USED* 


by 


Martin  Diamond 


September  7,  1971 


Department  of  Computer  Science 
University  of  Illinois  at  Urbana-Champaign 
Urbana,  Illinois  61801 


*This  work  was  supported  in  part  by  the  National  Science  Foundation  under 
Grant  No.  GJ812. 


i.  iwrroDucTioN 

The  solution  of  a  system  of  equations  AX  =  q  by  factorization 
techniques  can  be  separated  into  two  segments.   The  first  step  is  the 
definition  of  an  auxiliary  matrix,  B,  which  is  chosen  so  that  A  +  B  =  LU 
where  L  and  U  are  sparse  lower  and  upper  triangular  matrices  respectively. 
In  practical  computer  programs  the  elements  of  B  are  not  computed;  instead 
the  elements  of  L  and  U  are  computed.   Then  by  using  the  factors  L  and  U, 
systems  of  equations  whose  matrix  of  coefficients  is  A+B  can  be  easily 
solved. 

The  second  part  of  the  method  consists  of  an  algorithm  which 
uses  an  iteration  of  the  form 


(A+B)Xn+1  =  (A+B)Xn  -  (AXn-q),    n  =  0,  1,  2,  ... 


to  generate  a  sequence  of  vectors  {X  }  which  converges  to  the  solution  AX  =  q. 

The  convergence  of  the  iteration  is  dependent  upon  the  definition  of  the 
auxiliary  matrix  (through  the  definition  of  L  and  U)  and  a  sequence  of 
parameters. 

One  of  two  types  of  parameters  is  employed  in  the  iteration.   The 
first  type  multiplies  the  term  (AX  -q)  in  the  iteration  as  follows, 


(A+B)X  ,_  =  (A+B)X  -  t  (AX  -q).  (l.l) 

n+1   N   '   n    n   n 


an  analysis  of  this  parameter  and  an  algorithm  based  on  the  iteration  (l.l) 
can  be  found  in  Diamond  [  1] . 

The  second  set  of  parameters  may  be  used  to  alter  the  auxiliary 
:  matrix  for  each  n,  producing  the  iteration 

(A+B  )X  .  =  (A+B  )X   -  (AX  -q).  (1.2) 

n  n+1      n  n      n 

(Although  there  are  algorithms  which  use  a  sequence  of  auxiliary  matrices; 
jthere  has  been  no  rigorous  mathematical  basis  for  the  definition  of  an 

jiuxiliary  matrix,  or  a  sequence  of  matrices,  which  will  make  (1.2)  rapidly 

i. 

j  convergent. 


Stone  has  developed  in  [3]  an  algorithm  of  this  type  which  has 

proved  to  be  rapidly  convergent  on  a  number  of  test  problems.   The 

fundamental  idea  of  Stone's  factorizations  is  to  design  the  auxiliary 

matrix  so  that  the  terms  in  each  component  of  B  X  cancel.  When  B  X  and 

n  n  n 

B  X    are  small  it  is  intuitively  plausible  that  one  iteration  of  (1.2) 

will  yield  a  vector  which  is  approximately  equal  to  the  solution  of  AX  =  q. 

This  cancellation  is  the  basis  of  the  factorization;  however  it 
doesn't  explain  the  convergence  of  the  iteration.   Indeed,  experimental 
results  have  shown  that  repeatedly  using  the  auxiliary  matrix  B  _, 

(see  3*2  and  3*3  for  precise  definition)  which  produces  the  greatest 
canceling  yields  a  divergent  iteration. 

In  this  paper  we  explain  the  divergence  of  iteration  (1.2)  if 
the  auxiliary  matrix  B  ,  is  used  on  every  iteration.  We  then  present 

some  experimental  results  that  indicate  the  eigenvalue  of  the  iteration 
matrix  which  has  largest  modulus  grows  without  bound  as  the  dimension  of 
A  increases. 


3 

2.   ORIGIN  OF  THE  PROBLEM 


Consider  the  Dirchlet  problem 


-  st:  <ai»  IB  -  i:  ^  1-)  =  .w,  xed, 


u(X)  =  i|r(X),  for  XeBD  (2.1) 

where  x  =  (x,,  x  ),  D  is  the  interior  of  a  compact  region  with  "boundary 

dD  and  where  the  differential  operation  in  (2.1)  is  elliptic,  that  is 
a  (x)  and  a  (x)  are  strictly  positive  on  D.  A  common  method  of  solving 

(2.1)  is  to  approximate  the  differential  operator  with  a  difference 
operator.  A  system  of  linear  equations  results  and  the  problem  is  then 
to  find  a  solution  to  the  system  of  equations. 

The  system  of  equations  is  derived  as  follows.  Assume  D  lies 
in  the  first  quadrant;  cover  D  with  two  families  of  parallel  lines, 


x1  =  jh         j  =  1,  . . . ,  n1 


and 


x  =  kh         k  =  1,  . . . ,  n  . 


2 


Denote  the  points  of  intersection  (jh,  kh)  by  (j,k).   The  points  (j,k) 
are  called  the  mesh  points  or  gridpoints.   Let  dD,  be  the  polygon  formed 

by  joining  those  mesh  points  (i,  j)  for  which  one  of  (j"+l,k)  or  (j,k+l) 
is  not  in  D.   Similarly  let  D,  be  the  open  and  connected  set  bounded  by 

o\D,  .  At  each  gridpoint  of  D,  the  differential  operator  is  approximated  by 

a  difference  operator  and  the  equations  are  in  one-to-one  correspondence 

with  the  gridpoints. 

For  simplicity  of  notation  we  shall  consider  the  region  D  to  be 

the  unit  square  buonded  by  the  x,  and  x  axes  and  the  lines  x,  =  1  and 

T 
x0  =  1.   Letting  x  =  (xn  _.  x,  _,  ...  x  _,  x0  n,  ....  x   )  where 
2  °       1,1'   1,2'      n,  1'   2,1'    '   n,n 

n,  =  n  =  n,  we  obtain  the  matrix  equation 


AX  =  q  (2.2) 

defined  by  (AX)  .  ,  =  q.  .,  where  the  subscript  corresponds  to  the  gridpoint 

(j,k)  and  q.  ,  is  derived  from  q(jh,  kh)  and  the  boundary  conditions. 

The  difference  operator  can  be  chosen  so  that  the  matrix  A  is 
diagonally  dominant  and  positive  definite  with  a.  .  >  0  and  a.  .  <  0  for 

i  =  j.   One  method  for  deriving  such  an  A  is  to  approximate  ^ —  (a.  (x)  ^ — u( 

1        1 

by 


-2       h  h 

h  {ai(x+  -ei)[u(x)  -  u(x+hei)]  +  ai(x-  g  ei)  [u(x)  -  u(x-hei 


where  e-.  and  ep  are  unit  vectors  along  the  x,  and  x  axes  respectively.  Thi; 
defines  A  by 

(AX)  .  .  =B.  .X.  ,  -  +  D.  .  X.  _  .  +E.nx.n  +D.n  .x... 
'3,k    j,k  J,k-1    j,k  o-l, k    j,k  j,k    j+l,k  j+l,k 


+  Bj,k+lXj,k-KL'    0,k  =  1,  ...,  n,  where 


5j  k  =  "a2^h'  (k-2^11^ 


(2.3) 


\k  =  '^(ti^M, 


and 


Ej  k  =  a2(jh,  (k-|)h)  +  a2(jh,(k+|)h)  +  ai((j-|)h,kh) 


+  a1((j+2-)h,kh). 


When  j  =  1  or  k  =  1  we  define  D.  ,  =  0  or  B .  ,  =0  respectively. 

J,k        j,k 

The  following  properties  of  A  and  the  quantities  defined  in  (2.3)  will  be 

used  below. 


and 


(AX)  .  .    =  E  .  .  X .  .    +  B .  .  X .  .     ,+D.nX.nn     +  D  .    .    .  x .    _    . 


j,k+l  j,k+l 


B.   ,    <  0        ,         B.    .    =   0 
3,^  ~  j,l 


V'*°  >   v=o 


(2A) 


E.  .    >  -2(B.  .    +  D.  .  )  with  equality  holding  if  B.  .    ^  0 


and  D.  ,    /  0  for  j  =  1,   2,    ...,   n  and  k  =  1,    2,    . ..,   n. 


The  subscript  (j,k)  refers  to  the  gridpoint  (jh, kh). 


3.   THE  FACTORIZATIONS  OF  STONE 

Stone  has  defined  two  factorizations  both  of  which  depend  on  a 
parameter  a.   One  yields  a  symmetric  auxiliary  matrix  and  the  other  yields 
a  non symmetric  auxiliary  matrix. 

For  the  symmetric  case,  [4]  the  factors  L  and  U  of  A  +  B  are 
defined  as 


where 


and 


(LX)  .  .  =  b.  .  x.  ,  -,+cx     +  d.  X.  , 
3f*         3,k   D,k-1    g,k  o-l, k    j,k  j,k' 


(UK)  .  .  =  x.  ,  +  e.  ,x.  _  .  +  f .  .x.  .  _, 
;0,k    j,k    j,k  j+l,k    o,k  3,3*4-1' 


(3-D 


0,k  "   j,k"acj,k-l  j-l,k-l' 


c .   =  D.  -ab .  -  ,  e.     ,  (3-2) 

j,k   o,k   a-i,k  o-i,k-i' 


0,k    o,k  o,k-l    o,k  j-l,k    g,k    j.,k-l  j-l,k-l 


+  ab  .  i  ,  e  .  ,  ,  -, 
0-1,  k  o-l,  k-1, 


d .  .  e  .  ,  =  D . , _  .  -  ab  .  ,  e  .  .  , , 

0,k  j,k    J+l,k     j,k  3,k-l' 


df    =  B.  ,  ,,  -  ac  ,f.  t  ,  i 
0,k  o,k    o,k+l     g,k  o-l,k 


In  addition  we   define 


C .  .    =  Id  .  .  e .  ,    .    and  G .  ,    =  c .  ,  f .   .   ,  . 
3,k         o,k  j,k-l  j,k         oA  j-l,k 


7 


The  nonsymmetric  case  [3]  is  more  general  in  that  it  is  not 
necessary  to  have  A  defined  as  in  (2.^).   To  avoid  any  confusion  the 
original  matrix  will  "be  called  M  and  the  auxiliary  matrix  N.   M  may  be 
nonsyrametric  as  is  the  matrix  which  naturally  arises  in  the  solution 
of  a  problem  with  Neumann  boundary  conditions.  Let  M  be  defined  as 
follows : 


(MX).  ,  *=  B.  ,X.  ,  .  +  D.  '  X.  _  .  +  E.  .X.  . 
3,k    j,k  j,k_1    3,^   D-l,k    j,k  j,k 


,  j,k  j+i,k    j,k  3,k+l 

The  auxiliary  matrix  N  is  defined  by  factors  L  and  U  which  are  in  the  form 
(3.1)  just  as  in  the  symmetric  factorization.   Instead  of  (3.2)  however,  the 
elements  of  L  and  U  are  defined  by 

j,k  ""  j,k  "   j,k' 


cj,k  "   j,k  "   j,k' 


di  k  +  b-i  kfi  k-1  +  c-i  kei-l  k  =  Ei  k  +  ^i  k  +  ^-i  k> 


d.  ,e.  .  =  F.  .  -  qC  .  .  , 
3,k  3,3s    3,k     3,k' 


(3-3) 


3,k  3,k  .   3,k  "   3,k> 

where        C .  ,  =  b .  ,  e .  „  -  and  G .  n  =  c .  ,  f .  _  ,  . 
3,k    3,k  3, k-1      3,k    3,k  3-l>k 

Since  there  will  usually  be  no  confusion  as  to  whether  the  auxiliary 
matrix  is  defined  by  (3*2)  or  (3-3),  B  will  be  used  to  denote  either  of  them. 
With  L  and  U  of  the  form  (3-1)  A  +  B  (or  M  +  N)  is  given  by 


((A+B)X)  -    C(LU)X)  .>k  -  bJf^jM  +  \kej,k-l*j+l,k-l 


cj,kxj-l;k  j,k         j,k  d,k-l       Cj,kej-l,k'xj,k 


+  d.  , e .  .  x.  ,   .    +  c .  . f .   _   ,  x.   _   .    _    +  d.  .  f .  .  x.  . 

j,k  j,k  j+l,k         j,k  j-l,k  j-l,k+l         j,k  j,k  j,k  +  1 


k.      WHY  STONE'S  FACTORIZATIONS  ARE  EXPECTED  TO  CONVERGE 

In  this  section  we  explain  why  the  factorizations  defined  in 
(3*2)  and  (3*3)  with  a  =  1  are  expected  to  make  the  iteration 


(A+B)Xn+1  =  (A+B)Xn  -  (AXn-q) 


rapidly  convergent. 

The  matrix  M  which  multiplies  the  error  vector  E  =  X  -  X  on 

n       n 

each  iteration  is  called  the  iteration  matrix.   It  is  well  known  that  an 

iteration  is  convergent  if  the  spectral  radius  of  the  iteration  matrix  is 

less  than  one.  Moreover  as  the  spectral  radius  decreases,  the  rate  of 

convergence  increases.   Therefore  we  will  show  how  the  factorizations  (3*2) 

and  (3*3)  are  designed  to  produce  an  iteration  matrix  with  a  small  spectral 

radius . 

Throughout  the  following  the  eigenvalues  of  an  n  x  n  matrix  M 

will  he  denoted  by  \. (M)  with 


WM)  =  VM) ■>  VM)  •••  ^VM)  =  WM)- 


The  iteration  matrix  of 


(A+B)Xn+1  =  (A+B)Xn  -  (AXn-q)  (k.l) 


is  M  =  I  -  (A+B)  A.   The  spectral  radius  of  this  matrix  is  small  if  the 
\.  ((A+B)  A)  are  close  to  1.   Thus  the  iteration  is  rapidly  convergent  if 

(A+B)  A  is  approximately  equal  to  the  identity  I.   That  is  if 


lX-(A+B)  AX||  ig  small  fQr  all  x  ^  Q>  ([)_>2) 


PI 
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The  factorization  (3*3),  suggested  by  Stone,  attempts  to  satisfy 
(^.2)  by  making  ||(A+B)X  -  AX||  small  in  the  sense  that 


(BX).   =  o(h2)  (k.3) 


where  h  is  the  meshsize  of  the  grid  used  in  the  discretization  of  the 
differential  problem,  i.e.,  deriving  (2.2)  from  (2.1). 

To  see  that  BX  satisfies  the  condition  (k-3)   we  have  the  following. 

Let  x  =  (x_  ..,  ....  x   )x  where  the  x.  .  are  values  of  a  smooth  function 
1,1'    '   n,  n  j,k 

taken  at  the  gridpoints  (jh,kh).   If  B  is  the  matrix  given  by  (3«3)  and  C_ 


3,k 

definition  of  A  in  (2.1 

to  the  definition  of  A+B  in  (3-3)  and  (3-^)  yields 


and  G.   are  also  as  in  (3.3),  then  comparing  the  definition  of  A  in  (2.k) 


[(A+B)X]  .  .  =  (AX)  .  .  +  C.  .  [x.  in  ,  1    +   a(x.  .  -x.^-  ,  -x.  ,   ' 
J0,k      yj,k    j,k  L  j+l,k-l     j,k  j+l,k  j,k-l' 

+  G.n  [ x .  _  .  , _  +  a(x .  ,  -x .  ,  .  -x .   .)]. 
j,k   j-l,k+l     3,k  j-l,k  o,k+l/J 

V/hen  a  =  1,  Taylor's  Theorem  implies  the  bracketed  terms  of  (k.k)   satisfy 


2 

x       +  (x.  -x...  ,  -x.   _ )  =  0(h  ) 

j+l,k-l     j,k  j+l,k  j,k-l7 


and 


x.  n  .  .  +  (x.  ,  -x.  _  ,  -x.  ,  -, )  =  0(h  ). 

,]-l,k+l     j,k  g-l,k  j,k+l; 

Stone's  symmetric  factorization  (3.2)  satisfies  (^.2)  in  the  sense 
BX)  .   =  0(h).   This  also  follows  from  applying  Taylor's  Theorem  to 

(BX)    =  ((A+B)X-AX)    . 

Of*-  d)^~ 
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5.   WHY  STONE'S  FACTORIZATIONS  DIVERGE  WHEN  NO  PARAMETERS  ARE  USED 

The  meaning  of  small  used  above  is  "weak  since  (^.2)  and  (k.3) 
are  not  equivalent.   In  fact,  empirical  results  show  that  both  factorizations 
(3.2)  and  (3.3)  yield  a  B  such  that  p(l-(A+B)  A)  >  1  when  A  is  derived  from 
the  Laplacian  [a  (x)  =  a  (x)  =  1  in  (l.l)]  and  n  =  30.   Thus  even  though 

(BX)  .  .  =  0(h)  or  0(h2),  ||X-(A+B)-1AX||  is  not  small.   It  isn't  even  less  than 

1.  In  this  section  we  will  see  how  the  above  situation  is  possible. 

We  shall  begin  by  deriving  some  properties  of  the  eigenvalues  of 

(A+B)"^. 


Lemma  5.1: 

If  A  is  given  by  (2.k)   and  B  is  defined  by  either  of  Stone's 
factorizations,  (3*2)  or  (3-3)>  then  one  is  an  eigenvalues  of  (A+B)  A. 


Proof: 

If  B  is  given  by  (3»2)  or  (3*3)  then  the  first  row  and  column  are 
null.  Hence  B  has  a  zero  eigenvalue.   Let  w  be  the  corresponding  eigenvector. 

Then  (A+B)w  =  Aw  and  consequently  (A+B)  Aw  =  w  which  completes  the  proof. 


Theorem  ^.1: 

The  vector  w  is  an  eigenvector  of  (A+B)  A  iff  Bw  =  0  or  Aw  =  o"Bw 
where  <j  is  a  scalar  other  than  0  or  -1.  Furthermore  the  corresponding 

eigenvalue  is  -2=r  if  Bw  J   0  and  is  1  if  Bw  =  0.  (5-1) 

CT+-J- 


Proof : 

Suppose  Aw  =   CTBw       a  ^   -1,    0.      Then 


w  +  A     Bw=w  +  -w=    (- — )w. 
a  a   ' 
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Hence 


(A+B)_1Aw  =    (A+B)"1A(-^T)(w+A"1Bw) 

a+-L 


=    (-2T)(A+B)"1(Aw+Bw)   =  -2-w. 
CT+-L  a+J- 


Now  suppose    (A+B)     Aw  =  —rr  v,    a  4  ~lj    0.      Then 

cr+J- 


Aw  =  -MA+B)™ 

CT+1 


and  therefore  Aw  =  aBw. 

The  case  Bw  =  0  was  taken  care  of  in  Lemma  5.1  and  the  proof  is 
therefore  complete. 

For  simplicity  we  shall  now  consider  only  the  symmetric 
factorization  (3.2).   The  eigenvalues  \.    =  A..  (l-(A+B)~  A)  are  real  and 

have  magnitude  greater  than  one  iff 


°i  1 

A..  I  =  ll  -  =-  I  where  -2  <  a.    <  -—. 

i    '    cr.  +  1  '         —i—2 


If  A  is  defined  "by  (2.k)   and  (A+B)  is  defined  by  (3 .2)  then  it  is  easily 
shown  that  A..  ( (A+B)"^)  >  0,  see  Diamond  [1].   This  and  Theorem  5.1  imply 

ai  1 
T-T  >  0.   Thus  if  a-   <  -o  then  a.    +   1  must  also  be  negative  and  cr.  <  "!• 

CT .  +  J.  1      <d  1  1 

Hence  the  iteration  (^.1)  is  divergent  if  there  exists  some  w  such  that 

Aw  =  crBw,  -2  <C  cr  <;  -1.   Moreover  as  0   approaches  -1,  (l  -  — ~t)   approaches 
—   —  a  +  ± 

infinity  and  the  iteration  diverges  rapidly,  i.e.  the  error  grows  rapidly. 
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The  iteration  is  rapidly  convergent  "when  \\.  |  is  small  or,  in 
Dther  words,  when  \a-\    is  large.   The  fact  that  (Bw)  .   is  small  is  therefore 
Dnly  significant  if  it  is  small  relative  to  (Aw)    .   If  for  some  w  ||Aw||  is  small, 

say  0(h  ),  then  selecting  a  B  which  satisfies  (^.3)  does  not  guarantee  (^.2). 

rhis  difficulty  arises  on  even  the  simplest  problems  as  seen  from  the  following 

considerations.   It  is  well  known  (see  for  example  Frankel  [2]  that  for  A  the 

2    2 
(n  x  n  )  matrix  derived  from  the  Laplacian, 


X    .  (A)  = — 5  =   0(h2)  where  h  = 


min      2(n+l)2  '      '     '     "  n  +  1 


Thus  in  this  case  the  fact  that  B  as  defined  by  Stone  satisfies  (EX)    =  0(h~ 
or  0(h  )  doesn't  imply  that  ||x  -  (A+B)"1AX||  is  small  for  all  X  ^  0. 


Ik 


6.   EMPIRICAL  RESULTS  CONCERNING  X^^ 


In  this  section  we  shall  present  some  empirical  results  which 

indicate  \>r„,r((A+B)  A)  grows  without  "bound  as  the  dimension  of  A  increases, 
MAX 

n 
This  is  done  by  showing  the  existence  of  vectors w  ,  where  the  n  corresponds 

to  the  dimension  of  A,  such  that 


<  Aw  ,w  > 
lim  l =  co 

n-*»  <  (A+B)w  *,w 


Then  using  the  equality 


we  conclude    lim  ^MAY  ((A+B)  A)  =  0. 
n-*» 


n 
The  results  of  Section  5  suggest  the  vectors  w  will  satisfy 


Awn  =  a   Bw11  (6.2) 

n 


where 


lim  a     =  -1. 

n 
n-*» 


The  vectors  do  not  satisfy  (6.2).   However,  as  n  increases 


(Aw  )    S  -(Bwn)  ,  for  j  /  l,n  and  k  J-   l,n: 


and 
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<  Aw  ,    w     >  =  a     <  Bw  ,w     > 
'  n  ' 


•where  (6.3) 


lim  a     =   -1« 
n 

n-*» 


From  (6.3)  we  see  that 


.  .   n     n  ,. 
<  Aw  ,w    > 

1 

ffn 

<  (A+B)wn,wn  > 

CTn 

1  +  *n 

and  this  approaches  infinity  as  a  ^  -!•  Furthermore  the  vectors  w  do 
not  maximize 


<  AX,X  > 
<  (A+B)X,X  >  ' 


Instead  they  give  a  lower  bound  for 


$  <  (A^')X,X  > 

and  thus,  from  (6.1),  they  give  a  lower  hound  for  A.  .  ((A+B)  A). 

The  matrix  A  is  derived  from  the  Laplacian  and  Bis  the  auxiliary 
matrix  defined  in  (3*2)  with  a  =  1.   (The  results  also  apply  to  the  auxiliary 
matrix  defined  in  (3*3)  since  the  elements  of  the  two  matrices  approach  the 
same  values  as  the  order  of  A  increases.   This  can  be  seen  by  noting  that  the 
proof  of  Theorem  5-1  in  [1]  can  be  applied  to  the  elements  of  the  nonsymmetric 
factorization  as  well  as  the  symmetric  factorization. ) 

A  number  of  different  vectors  were  tried  on  a  trial  and  error  basis 

in  attempting  to  find  a  w  such  that  <  Aw  ,    w  >  =  -  <  Bw  ,  w  >.   The  process 
of  selecting  vectors  was  started  by  using  the  algorithm  of  Diamond  [1]  to 
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calculate  the  eigenvector  w  corresponding  to  the  largest  eigenvalue  of 

(A+B)  A.   Next  a  function  f(x  ,  x  )  was  found  such  that  f(jh,  kh)  =  w   . 

A  vector  F  =  (f    ,  ...,    f    ),   where  f     =  f(jh,kh)  and  h  =  l/(l+n), 

was  defined  and  used  as  the  first  trial  vector.  A  more  or  less  random 

procedure  was  employed  to  alter  the  trial  vector  F.  and  form  the  succeeding 

n 
vector  F.  _ •   The  final  vector  is  denoted  by  w  ,    where  as  before  the  n 

corresponds  to  the  dimension  of  A. 

From  the  eigenvector,  which  was  calculated  to  select  the  first 
trial  vector,  the  following  information  was  inferred  about  the  function 
f(xx,x2). 


The  function  should  have  a  maximum  at  the  center  of  the 

rl 
,2i 


,1   In 

region,  (-,  -J. 


2.  The  function  is  symmetric  with  respect  to  the  lines  x  =  -x  , 

1       i        .1 

X_  -  r>)      Xp  -  1   anC!-  X"I   _  O* 

3.  The  function  approaches  zero  as  (x  ,  xp)  approaches  the 

boundary. 
k.      The  function  decreases  faster  along  the  line  x  =  x  than 

along  x  =  1  -  x  . 


Functions  with  the  above  properties  were  tried  with  the  additional 
possibility  that  they  were  periodic  in  the  x  and  x  directions  with  period 

Tr/m  ,   m  an  integer.   In  this  case  the  function  exhibited  the  above  propert: 

in  each  period. 

n   n  n   n 

The  intention  was  to  determine  <  Aw  ,  w  >  and  <  Bw  ,  w  >  when  n 

was  very  large.   That  is  when  B  may  be  partitioned  as 


B  = 


Bl,l   Bl, 


'2,  1   \,  2/ 


IT 


where  B„  0   is  an  {W   X  HT)  submatrix  with  almost  constant  diagonals.   Using 
the  limiting  values  of  the  elements  of  B  as  derived  in  [ 1]  we  see  that 


B2,2  "  B2,2 


Instead  of  computing  <  Bw  ,w  >  the  quantity  <  Bp  ,  pw  ,w  >  was 
computed.  We  have  <  B  ,      w  ,   w  >  =  <  B  ,   w  ,  w  >  =  <  Bw  ,    w  >  where 


2   '    2 


2'    2 


n 
w     = 


w 


and  0  is   an   (n  -IT) -component  zero  vector „      Similarly  <  Awn,wn  >  was  computed  from 


.An       n  ^        ^  .  N       N         , 

<  Aw  ,   w     >  =  <  A     _w  _,   w     >  where 


A 


\l       Al,2 


k2, 1       A2, 2 


Then  instead  of  computing 


an  : 

/      A       n 

<  Aw  , 

<  Bw  , 

n  ^ 
w     > 

n 
w     > 

<  A2  gw  ,    w     > 

<  B2  2w  ,    w     > 

the  quantity 
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^       <  A2  2w  ,  w  > 

°N  =  <  w,  ^  > 


was  computed. 


N 
The  vector  w  ,  as  described  above,  was  chosen  to  be 


v  kN  =  exp  {-(|jh-||  +  |kh-|l  +  ±f   [2  -  cos( (jh-kh)^|)]} 

cos2°((jh-kh)^|). 

-x         -x- 
For  N  =  100  and  N  =  250  the  values  of  a       were  a_    =  -1.0918  and 

-X  -x- 

a„.-~  =  -1.066l.   Values  of  N  between  100  and  250  produced  values  of  u_T 
2pU  N 

-X  "X* 

which  varied  monotonicly  from  cr-,m  to  Open  •   This  indicates  that  as  N 

-x 
increases  cr   gets  close  to  -1  and 


1  +  a  -,     * 

n    1  +  aN 

grows  without  bound.      Thus  X   .    ((A+B)     A)   grows  without  bound  as  the 
dimension  of  A  increases. 
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